#include <sfv/array1d.h>

#include <stack>
#include <deque>

namespace sfv {

template< typename MatrixType, typename Array1, typename Array2 >
int detect_clusters(const MatrixType& graph, Array1& membership,
		Array2& csize, size_t *no) {

	typedef typename MatrixType::index_type IndexType;
	typedef typename MatrixType::value_type ValueType;

	long int no_of_nodes = graph.num_rows;
	long int act_cluster_size = 0, no_of_clusters = 1;

	std::deque<long int> q;

	sfv::array1d<char> already_added(no_of_nodes,0);

	/* Memory for result, csize is dynamically allocated */
	membership.resize(no_of_nodes);
	csize.clear();

	/* The algorithm */

	for (long int first_node = 0; first_node < no_of_nodes; ++first_node) {
		if (already_added[first_node] == 1)
			continue;

		already_added[first_node] = 1;
		act_cluster_size = 1;
		membership[first_node] = no_of_clusters - 1;
		q.push_back(first_node);

		while (!q.empty()) {
			long int actnode = q.front();
			q.pop_front();

			size_t neino = graph.row_offsets[actnode+1] - graph.row_offsets[actnode];
			sfv::array1d<IndexType> neis(neino);
			IndexType base = graph.row_offsets[actnode];
			for( size_t i = 0; i < neino; i++ )
				neis[i] = graph.column_indices[base+i];

			for ( size_t i = 0; i < neis.size(); i++) {
				long int neighbor = neis[i];
				if (already_added[neighbor] == 1) {
					continue;
				}
				q.push_back(neighbor);
				already_added[neighbor] = 1;
				act_cluster_size++;
				membership[neighbor] = no_of_clusters - 1;
			}
		}
		no_of_clusters++;
		csize.push_back(act_cluster_size);
	}

	/* Cleaning up */

	*no = no_of_clusters - 1;

	return 0;
}

} // end namespace sfv
